# 加载需要的包
library(DESeq2)
library(tidyverse)

# 读取文件
patient_info <- read.table("D:\\RLocal\\Class7\\RNAseq_Data\\clin_inf.csv", row.names = 1, header = TRUE)
expr_data <- read.table("D:\\RLocal\\Class7\\RNAseq_Data\\count.csv", row.names = 1, header = TRUE)

# 先把列名统一
col_names <- colnames(expr_data)
# 对 PTB 开头的列名，去除最后一位
col_names <- ifelse(grepl("^PTB", col_names),
                    substring(col_names, 1, nchar(col_names) - 1),
                    col_names)
# 对 X 开头的列名，去掉最后一位，并把 X 换成 XYA
col_names <- ifelse(grepl("^X", col_names),
                    paste0("XYA", substring(col_names, 2, nchar(col_names) - 1)),
                    col_names)
# 应用修改后的列名
colnames(expr_data) <- col_names


# 创建分组变量
patient_info$AgeGroup <- ifelse(patient_info$年龄 >= 65, "65up", "under65")

# 保证表达矩阵与分组信息对齐
common_samples <- intersect(colnames(expr_data), rownames(patient_info))

expr_data <- expr_data[, common_samples]
patient_info <- patient_info[common_samples, ]

# 四舍五入转换为整数
expr_data_int <- round(expr_data)

# 确保是整数矩阵（因为round后结果是numeric）
expr_data_int <- apply(expr_data_int, 2, as.integer)

# 重新变成data.frame，保留行名
expr_data_int <- as.data.frame(expr_data_int)
rownames(expr_data_int) <- rownames(expr_data)

# DESeq2分析
dds <- DESeqDataSetFromMatrix(countData = expr_data_int,
                              colData = patient_info,
                              design = ~ AgeGroup)
# 过滤低表达基因
keep <- rowSums(counts(dds)) > 10
dds <- dds[keep, ]

# 差异表达分析
dds <- DESeq(dds)
res <- results(dds, contrast = c("AgeGroup", "65up", "under65"))

# 保存结果
write.csv(as.data.frame(res), "DEG_by_AgeGroup.csv")